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Abstract 

This report contains a detailed theoretical de- 
scription of an all-purpose, interactive curve-fitting 
routine that is based on P. R. Bevington’s descrip- 
tion of the quadratic expansion of the x 2 statistic. 
The method is implemented in the associated inter- 
active, graphics-based computer program. 

The Taylor’s expansion of x 2 is first introduced, 
and justifications for retaining only the first term 
are presented. From the expansion, a set of n 
simultaneous linear equations are derived, which are 
solved by matrix algebra. 

A brief description of the code is presented along 
with a limited number of changes that are required to 
customize the program for a particular task. To eval- 
uate the performance of the method and the good- 
ness of nonlinear curve fitting, two typical engineer- 
ing problems are examined and the graphical output 
and the tabular output of each are discussed. A com- 
plete listing of the entire package is included as an 
appendix. 


Symbols 

a 3 

Sdjy Afly 


y(xi),y(x) 

a jk 


0k 

£ jk 




Xo 

xl 


Subscripts: 

i 


3 

k 


v 


coefficient of fitting function 

coefficient of continuous and dis- 
crete differential of ay, respectively 

experimental data 

fitting function 

algebraic notation for symmetric 
matrix 

algebraic notation for row matrix 

inverse matrix of ay*. 

uncertainty in data 

global chi-square 

first term in the expansion of x 2 

reduced chi-square 

index of experimental data 

index of coefficient of fitting func- 
tion, also row index of a symmetric 
matrix 

column index of a symmetric matrix 
number of degrees of freedom 


Introduction 

In any area of engineering or physical science, 
suggested analytical models are accepted only when 


good statistical correlation exists with a set of exper- 
imentally measured values. The correlation is often 
measured by fitting the mathematical model to a set 
of experimental data. 

Two common methods for fitting data are mov- 
ing averages and least-squares fit. In the moving av- 
erages method, each data point is replaced by the 
average of itself and n neighboring points on either 
side of it. The advantage of this method is that it 
is rather easy to program. One disadvantage is un- 
equal smoothing of the first and the last data points 
compared with the rest of the data set because of the 
lack of neighbors on both sides. Another, more im- 
portant, disadvantage is that the smoothing process 
is strictly an averaging one and does not produce any 
analytical representation of the smoothed data. 

In the least-squares method, a user-specified fit- 
ting function is utilized in such a way to minimize 
the sum of the squares of distances between the data 
points and the fitting curve. The advantages of this 
method are that it permits the generation of statis- 
tical information on the goodness of the fit and does 
not require the data to be collected at regular in- 
tervals. The disadvantages are that the method as- 
sumes that the basic form of the smoothing equation 
is known and also, since it is a global procedure, it 
may be disproportionately biased by a few bad data 
points, which will twist the shape of the fit to spread 
the error over the entire data set. 

Considering the advantages of the least-squares 
fitting method and the decreasing expense of com- 
putation time, it is often desirable to have a consol- 
idated software package in the form of a single com- 
puter program to perform nonlinear curve fitting to a 
given set of data. This approach should provide the 
user with statistical information such as goodness of 
fit and estimated values of parameters that produce 
the highest degree of correlation between the exper- 
imental data and the mathematical model. 

The purpose of this paper is to furnish such a soft- 
ware package. The section “Fitting Algorithm De- 
scription” describes the mathematical formulation of 
the quadratic expansion of x 2 > which fundamentally 
follows the work of Bevington (ref. 1) and in many 
cases closely parallels his discussion. The section 
“Program Description” briefly describes the modu- 
lar characteristics of the program and its associated 
subroutines and function subprograms. These pro- 
gram elements are formulated around a nonlinear 
optimization algorithm that calculates the best sta- 
tistically weighted values of the parameters of the 
fitting function and the x 2 that is to be minimized. 
The program needs as input the mathematical form 
of the fitting function and the initial values of the 
parameters to be estimated. The “Notes to Users” 


section describes the limited changes a user must 
make to implement the program for a particular ap- 
plication. The section “Sample Cases” describes two 
sample cases. 

Fitting Algorithm Description 

Consider the function y(x) with parameters ay. 
For example, y(x ) can be an exponentially decaying 
sinusoidal function, plus a constant, of the form 


the x 2 hypersurface and to use this function to locate 
the minimum. 

Description of x 2 Expansion 

Consider the linear terms of a Taylor expansion 
of x 2 as a function of parameters ay 

*w„+J;(^W) (5 ) 


y(x) = die a2X cos(dgx + a±) + 05 (la) 

or, a double Gaussian function, plus a quadratic, of 
the form 


y(x) = die 2 V “ 3 / 

_lf x ~ a 5 A 2 

+ d±e 2 V a 6 / + aj + dgx + agx 2 (lb) 

or some other function such that some of the param- 
eters cannot be separated into different terms of a 
sum. 

Bevington (ref. 1 ) defines x 2 > a measure of the 
goodness of the fit, as 


where 8dj are the increments in dj required to reach 
the point at which y{x) and x 2 are to be evaluated. 
The Xq is the starting value of x 2 at the point where 
the value of y(x) is yo(x) such that 

xo = E 1 ^ 2 1 y* ~ yo( x i)? J ( 6a ) 

and 

2/o(z) = y{x, aio, a 2 0,..., a n0 ) (6b) 

Since the optimum values for dj are defined 
through the minimization of x 2 with respect to ay, 
then 


X 2 = E|^I!'.-!'W] 2 } < 2 ) 

where a?, the uncertainties in the data points y t -, is 
defined as 

a i = “ E “ &) 2 ( 3 ) 

. - 
J = 1 

According to the method of least squares, the 
simultaneous minimization of x 2 with respect to each 
of the parameters produces the optimum values of 
parameters dj as 

Because of the difficulty in deriving an analytical 
expression to calculate the parameters of y(x), x 2 ls 
considered as a continuous function of n parameters 
dj describing a hypersurface in a space of n + 1 di- 
mensions, where ay, j = 1, 2,..., n, are the abscissa 
and x 2 is the ordinate. This space is searched to 
locate the minimum value of x 2 - 

In the present paper the search is accomplished 
through the expansion of x 2 by using an analytical 
expression for the variation of x 2 to map its variation 
with respect to parameters ay. The goal will be to 
find an approximate analytical function describing 


dx 2 _ dxo 

da^ dajt 



o2 2 
£ Xo 

daj dak 



= 0 


(k = 1, 2 ,..., n) 

(7) 


A set of n simultaneous linear equations in 6aj are 
obtained, which algebraically can be written as 


Pk = E ( Sa j a jk ) (* = 1, 2, n ) (8a) 

3 = 1 


where 


/j.=_I£xo a . - 1 

- 2 da k Jk ~ 2 daj da k 


(8b) 


One way of looking at equation (8b) is to state 
that x 2 through the first-order expansion is approx- 
imated by a parabolic surface. This is verified by a 
second-order Taylor expansion of x 2 as a function of 



which is a second-order function with respect to 8dj 
and describes a parabolic hypersurface. 
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Equation (9) indicates that the optimum values 
of 6a j for which x 2 is a minimum are obtained by 
requiring that the derivatives with respect to aj be 
zero. Thus, 



do k & a k 


+ 


n / ^2 2 

E l d Xo 

3 _ y \^ a j ^ a k 



= 0 


{k= 1 > 2 ,..., n) 
( 10 ) 


which is the same as equation (7). 

The method of quadratic expansion is accurate 
and precise if the minimum is close to the starting 
point in such a way that higher order terms in equa- 
tion (9) can be neglected. But, if the starting point 
is not close enough, the parabolic approximation of 
X 2 hypersurface is generally not valid, and in the di- 
rection of increasing 6a j the result will be in error. 
Hence to achieve convergence the algorithm requires 
meaningful initial estimates for aj. The initial esti- 
mates can often be obtained by visual inspection of 
data. 


Description of Computational Method 

The analytical methods of the previous section 
can be used for computational purposes by recogniz- 
ing that a matrix inversion operation will yield the 
solution of equation (8) as 


n 

$ a j = ^2 (Pk € jk) 


k= 1 


( 11 ) 


where = c*^ 1 , and the computation of equa- 
tion (8b) can be approximated by calculating the 
variation of x 2 near the starting point Xo an( l us- 
ing the standard finite difference equations of first, 
second, and cross product derivatives of Xo with re ” 
spect to daj and daj da ^ as 



2Aa • ^ + ^ a j' a k) 

- Xo (oj- AapOfc)] 


(12a) 


J j 

-2\o {a j,a k ) 

+ Xo(aj “ Aa i»°fc)] ( 12b ) 


<9^x 2 1 r 

da 3 da k * AoyAofe i X ° ^ j + Aaj ’ Uk + A ° fc) 

- Xo i a j + a k ) - xo (aj, a k + A a k ) 

+ Xo(«j’ a fc)] ( 12c ) 

Finally, the quantity v , the number of degrees of 
freedom left after fitting N data points to a function 
of n + 1 parameters, is defined as 

v = N-n-1 (13) 

Therefore, for u degrees of freedom, the quantity xt, 
the reduced chi-square, is defined as 



X^ will be used in the computations where N and n 
have specific numerical values. 

Program Description 

The program evolved from the idea of having an 
interactive package that requires minimum modifi- 
cation by the user. The main program and each 
subroutine or function subprogram begins with a 
description of its purpose and a definition of the 
variables used. The program is 882 lines long 
and is written in FORTRAN 77. It was devel- 
oped on the CDC® CYBER 750 scalar mainframe 
under the NOS 2.3 Level 617 Operating System 
and requires a minimum of 12400s 60-bit words 
of storage. The entire package is divided into a 
main program (NLNFIT), five subroutines (CHI- 
FIT, MATINV, PRETTY, CHAR, ERRBAR), and 
three function subprograms (FCHISQ, FUNCTN, 
TEXP), with the main program (NLNFIT) contain- 
ing all EXTERNAL Tektronix (PLOT- 10) CALLS 
(refs. 2 and 3) and Character Generator System 
(C.G.S.) CALLS (ref. 4). Subroutines CHIFIT and 
MATINV and function subprogram FCHISQ were 
originally developed in reference 1 and were modified 
by the authors. A brief description of the function of 
main module NLNFIT follows. 

Main Program (NLNFIT) 

NLNFIT assumes that the input data file named 
“RAWDAT” is written on logical unit 1 (LU = 1) as 
is specified by the PARAMETER statement. This 
can easily be changed to another suitable value if 
LU = 1 is a reserved unit. 

For the sake of transportability, the data file is 
limited to only four sets of input. The first card 
is an integer specifying the number of data pairs 
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and is optionally set at 200. The second card is an 
integer flag with values +1, 0, or -1, depending on 
whether the input data are to be weighted or not. For 
instrumental weight, where the uncertainty in each 
measurement of generally comes from fluctuations 
in repeated readings of instrumental scale, the input 
weight flag should be set to +1. The choice of 
instrumental weight requires that the user input data 
points (x^, yi) and uncertainty A y^. If it is decided 
not to weight the input data, integer flag 0 must be 
chosen. For statistical weight, where it is generally 
true that the uncertainty in each measurement yi 
is proportional to | y ^l” 1 and therefore the standard 
deviations ay associated with these measurements 
cannot be considered equal over any reasonable range 
of values, an integer flag of -1 must be chosen. The 
third card is the form of the fitting equation and will 
be read by the main module in an 80 A 1 format. The 
actual data pairs are the fourth input and are read 
in the form (x;, yi) for no weight or statistical weight 
or (x^, yi, Ayi) for instrumental weight. 

Program Execution 

The execution begins with the program asking 
for initial estimates of ay, j = 1, 2,... n, where n 
is the number of parameters. The output begins by 
informing the user if he has exceeded the limits of 
data pairs in the PARAMETER statement. If the 
limits have not been exceeded, the program displays 
the number of data pairs, the mathematical form 
of the fitting function, and the values of initial ay 
estimates. 

At this stage NLNFIT calls SUBROUTINE 
CHIFIT. This subroutine uses a quadratic expansion 
of the x 2 statistic to make a least-squares approxi- 
mation to the fitting function. 

During each iteration of CHIFIT (optionally set 
at 20 in the PARAMETER statement), NLNFIT dis- 
plays the iteration index and the value of the reduced 
CHISQR x 2 • The iteration continues until the differ- 
ence between two consecutive values of CHISQR is 
less than 1 percent or maximum iteration is achieved; 
in either case the final iteration index, values of ay, 
and averaged differences between yi and y(xi) are 
displayed, and the user is asked whether he wishes 
to see input values (x*,t/;) versus fitted y(xj). 

At this stage the user, if equipped with a Tek- 
tronix graphics terminal, is asked if he wishes to plot 
the input values of (x*,^) and the fitted t/(x*). If 
the answer is positive, a series of questions concern- 
ing the type of plot are asked. 


Notes to Users 

This section describes what changes a user must 
make to each routine (appendix A) to use the pro- 
gram for a different fitting function. 

NLNFIT 

The PARAMETER statement is the only change 
that is required for the main program. In the PA- 
RAMETER statement, II indicates the maximum 
number of data pairs, JJ must always be 4*11, KK 
is the maximum number of characters in the X and 
Y title statements, LL is the number of ay, IBAUD 
is the baud/ 10 rate of graphics display device, ITER 
is the maximum number of iterations allowed, and 
LU = 1 is the logical unit for input data. 

CHIFIT 

In CHIFIT, only the value of LL in the PARAM- 
ETER statement must be changed. 

FUNCTN 

In FUNCTN, the value of LL in the PARAM- 
ETER statement and the form of the FUNCTN 
statement must be changed. 

MATINY 

In MATINV, only the value of LL in the PARAM- 
ETER statement must be changed. 

Sample Cases 

Two sample cases in classical and fluid mechanics, 
weighted statistically (-1) and instrumentally (+1), 
respectively, are analyzed with the program package. 
Each case is described below, and its computer out- 
put is given as an appendix. 

Sample Case 1: Classical 

Mechanics — Physical Pendulum 

The circles in figure 1 are 166 data pairs obtained 
through an 8-bit A/D converter in a pendulum cali- 
bration test conducted by the authors. 

A 5-parameter nonlinear fitting function of the 
form 

A(t) ~ Aie -t /* ml cos(ujt + S) + (15) 

was applied to the data. Equation (15) is similar to 
equation (la), with a 2 = w and 6 as angular 
frequency and phase, and A 2 as contribution due to 
damping factors such as the frictional forces in the 
support bearings. The solid line is the best fit to the 
data. This particular functional form (eq. (15)), with 
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i 

I 

I 

I 

! initial dj estimates listed in appendix B, produced 
~ 0.12 in six iterations. Appendix B lists the 
interactive session for sample case 1. 


Sample Case 2: Fluid Mechanics — Far-Field 

Wind-Tunnel Pressure Analysis 

The circles in figure 2 are 22 data pairs repre- 
senting the nondimensional pressure coefficients mea- 
sured near the top wall of a two-dimensional wind 
tunnel. A 6-in-chord airfoil model was mounted 
on the tunnel centerline between x = — 3 in. and 
x = +3 in. The variation of the data is the result 
of the expansion of the flow about the model and 
a flow angularity probe inserted in the airstream at 
x = 6 in. near the top wall. The data were measured 
approximately 3.5 chord lengths above the model. 

A 9-parameter nonlinear fitting function of the 
form 


A(x) = Aie 2 V CT i / 


+ A±e 


(16) 


was applied to the data. Equation (16) is similar to 
equation (lb), with 02 = m, 03 = crj, 05 = 


and 06 = ^2 as the mean ft and standard deviation 
a of each Gaussian peak, and Aj, Ag, and .4 9 are 
the background contributions due to the undisturbed 
flow in the tunnel in the absence of the airfoil. The 
solid line is the best fit to the data. This particular 
functional form (eq. (16)), with initial aj estimates 
listed in appendix C, produced xi 0.69 in four 
iterations. Appendix C lists the interactive session 
for sample case 2 with initial data listed as X-DATA, 
Y-DATA, and fitted data listed as YFIT. 

Concluding Remarks 

The theoretical description of an all-purpose 
curve-fitting routine based on quadratic expansion 
of x 2 was presented. Taylor’s expansion of x 2 was 
introduced, and from the expansion a set of n simul- 
taneous linear equations were derived and solved by 
matrix algebra. The associated interactive, graphics- 
based computer program and sample cases indicated 
the relatively fast convergence rate of the method. 
Guidelines on how to customize the program for a 
particular task were given and fully described. 


NASA Langley Research Center 
Hampton, Virginia 23665-5225 
September 8, 1987 
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Appendix A 

Program Listing of Nonlinear Fitting Program NLNFIT 

Appendix A contains the program listing of the nonlinear fitting program NLNFIT, which consists of the 
main program NLNFIT, five subroutines CHIFIT, MATINV, PRETTY, CHAR, ERRBAR, and three function 
subprograms FCHISQ, FUNCTN, TEXP. 


7 



o o 


C PROGRAM NLNFIT 

C 

C PURPOSE 

C MAIN PROGRAM TO MAKE A LEAST SQUARES FIT TO A NON-LINEAR 

C FUNCTION WITH A QUADRATIC EXPANSION OF CHI. SQUARE 

C 

C DESCRIPTION OF PARAMETERS 

C II - MAX. NO. OF DATA POINTS <200) 

C JJ - 4 TIMES THE NUMBER OF DATA POINTS, USED FOR PLOTTING 

C A SMOOTH FIT THROUGH DATA POINTS 

C KK - MAX. NO. OF ALPHABETIC CHARACTERS IN TITLE STATEMENTS 

C LU - LOGIC UNIT OF I/O FOR INPUT DATA FILE 

C LL - NO. OF COEFFICIENTS OF FITTING FUNCTION 

C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE 

C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE 

C XDATA - DUMMY ARRAY TO STORE INDEPENDENT DATA POINTS 

C YDATA - DUMMY ARRAY TO STORE DEPENDENT DATA POINTS 

C YBU - DUMMY ARRAY TO STORE SIGMAY 

C YBD - DUMMY ARRAY TO STORE SIGMAY 

C YFIT- ARRAY OF CALCULATED VALUES OF Y 

C SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS 

C A - ARRAY OF PARAMETERS 

C SIGMAA - ARRAY OF STANDARD DEVIATIONS FOR PARAMETERS A 

C DELTAA - ARRAY OF INCREMENTS FOR PARAMETERS A 

C YLABEL - ARRAY TO STORE TITLE OF Y-AXIS 

C XLABEL - ARRAY TO STORE TITLE OF X-AXIS 

C IXLAB - DUMMY ARRAY FOR X-TITLE 

C IYLAB - DUMMY ARRAY FOR Y-TITLE 

C TITLE - ARRAY TO STORE FITTING FUNCTION 

C NPTS - NUMBER OF PAIRS OF DATA POINTS 

C MODE - DETERMINES METHOD OF WEIGHTING LEAST SQUARES FIT 

C +1 < INSTRUMENTAL > WEIGHT< I )«1. /SIGMAY< I >**2 

C O (NO. WEIGHT) WEIGHT ( I >«1. O 

C -1 (STATISTICAL) WEIGHT( I >-l. /Y< I > 

C NTS RMS - NUMBER OF PARAMETERS 

C CHISQR - REDUCED CHI. SQUARE FOR FIT 

C IBAUD - BAUD/10 RATE OF GRAPHICS DISPLAY DEVICE 

C ITER - NO. OF ITERATIONS TO CONVERGE (20) 

C 

PROGRAM NLNFIT ( INPUT, OUTPUT, TAPE5-INPUT, TAPE6-0UTPUT ) 

C 

PARAMETER ( 1 1-200, JJ-4*II, KK-BO, ITER-20, LU-1, LL-9, IBAUD-960) 


DIMENSION X(JJ) , Y(UJ), XDATA( JJ+1 ), YDATA(JJ+1>, YFIT ( JJ), SIGMAY <UJ) 
DIMENSION YBU< JJ+1 >, YBD (JJ+1 > 

DIMENSION A(LL) , SIGMAA(LL), DELTAA (LL) 

DIMENSION YLABEL(KK) , XLABEL (KK), IXLAB (KK) , I YLAB (KK ) , TITLE(KK) 

C 

WRITE (6, 10) 

C 

OPE N < UN I T-LU , AC CESS- ' SEQUENT I AL ' , F I LE= ' R AWD AT ' > 

C 

REWIND LU 
C 

READ(LU> 20) M»TS 
READ (LU, 20) MODE 
READ(LU, 270 ) <TI TLE< I ) , I-l.KK) 
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I 

I 

I 

I 

c 

DO 30 I = 1,NPTS 

IF ( (MODE. EQ. 0>. OR. (MODE. EQ. -1) > THEN 
READ(LU, *) X(I>.Y(I) 

ELSE 

IF (MODE. EQ. 1 ) THEN 
READ(LU.*) X(I>, Y(I),SIGMAY(I) 

END IF 
END IF 

I 30 CONTINUE 

C 

CLOSE (UNIT=LU) 

C 

IF(NPTS. GT. II) THEN 
! WRITE (6, 40) 

i STOP 

END IF 

WRITE(6, 50) NPTS 

NTERMS-LL 

KEEN=0 

j WR I TE (6, 60 ) 

WRITE (6, 270) (TITLE ( I )» 1=1, KK) 

WRITE (6, 70) NTERMS 

C 

DO 80 1 = 1, NTERMS 
WRITE(6, 90) I 
READ (5, *) A(I) 
i 80 CONTINUE 

! c 

^ WRI TE ( 6, 100) 

I WRITE(6. 110) (I, A(I)» 1 = 1, NTERMS) 

WRITE (6, 120) 

C 

DO 130 1=1, NTERMS 
DELTAA( I )=A( I)*. 01 
130 CONTINUE 
C 

K0UNT=0 

C 

C BEGIN ITERATION 

C 

DO 140 K-l. ITER 

CALL CHIFIT ( X, Y, SIGMAY, NPTS, NTERMS, MODE, A, DELTAA, SIGMAA, 
- YFIT. CHISQR) 

WRITE(6, 150) K, CHISQR 

IF ( K. GT. 1) THEN 

GO TO 160 

END IF 

SAVE=CHISQR 

K0UNT=1 

GO TO 140 

160 XCH I=CHISGR— SAVE 

IF ( ABS( XCHI ) . LT. 0. 01 ) THEN 
WRI TE(6, 175) 

GO TO 180 
END IF 
SAVE=CHISQR 
KOU NT =KOUNT + 1 
140 CONTINUE 
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180 KOUNT-KOUNT+1 

WRITE<6, 170) KOUNT 
WRITE <6. 190) 

WRITE <6, 270) (TITLE ( I ), 1=1, KK) 
WRITE (6< 45) 

C 

DO 200 1 = 1, N TER MS 
WRITE<6, 210) I, A(I) 

200 CONTINUE 
C 

WRITE<6, 220) CHISQR 
WRITE <6, 230) 

READ < 5. *) I ANSI 
IF< I ANSI. EQ. 1 ) THEN 
GO TO 240 
END IF 

490 WRITE<6, 250) 

READ< 5, *) IANS2 

IF( IANS2. EQ. 0) THEN 

STOP 

END IF 

WRITE<6, 260) 

READ<5, 270) <YL ABEL U ), 1 = 1 , KK > 
WRITE (6, 280) 

READ < 5, 270 ) <XL ABEL < I ) , 1 = 1 , KK ) 
WRITE (6, 290) 

READ < 5,*) IANS 
IF ( IANS. EQ. 2) THEN 
GO TO 300 
END IF 
GO TO 310 
300 WRITE (6, 320) 

READ( 5, *> INSL 
310 CONTINUE 

WR I TE ( 6, 330) 

READ<5, *) IANSRT 

IF( IANSRT. EQ. 0) THEN 

GO TO 340 

END IF 

WRITE (6, 350) 

READ< 5 , *) ISYMB 
340 CONTINUE 

WR I TE ( 6, 360) 

READ (5, *) NSETX 

IF ( NSETX. NE. 1 ) THEN 

GO TO 370 

END IF 

WR I TE ( 6, 380) 

READ< 5, *) XMIN, XMAX 
370 WRITE<6, 390) 

READ (5,*) NSETY 
IF<NSETY. NE. 1) THEN 
GO TO 400 
END IF 
WRITE (6, 410) 

READ (5,*) YMIN, YMAX 
400 CONTINUE 
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START OF TEKTRONIX PLOT-10 GRAPHICS CALLS 

CALL INITT < I BAUD) 

CALL BINITT 
CALL XNEAT(l) 

CALL YNEAT(l) 

XDATA< 1 ) “FLOAT ( 4*NPTS) 

YDATA( 1 ) “FLOAT < 4*NPTS) 

YBU ( 1 ) “FLOAT (4*NPTS ) 

YBD( 1 )=FL0AT(4*NPTS> 

FILL DUMMY ARRAY DATA POINTS 

DO 420 I“2» 4#NPTS+1» 4 
KEEN*KEEN+1 
XDATA( I )“X < KEEN ) 

XDATA( 1+1 )«X (KEEN) 

XDATA( I+2)“X (KEEN) 

XDATA( I+3)=X(KEEN) 

YDATA(I )“Y(KEEN) 

YDATA < I +1 >*Y (KEEN) 

YDATA( I+2)“Y(KEEN) 

YDATA< I+3)“Y(KEEN> 

IF (MODE. EG. 1 ) THEN 
YBD ( I )»Y(KEEN)-SIGMAY(KEEN) 

YBD ( 1+1 ) “Y ( KEEN ) -S I GMA Y ( KEEN ) 

YBD ( I +2 ) “Y ( KEEN ) -S I GMAY ( KEEN ) 

YBD ( I +3 ) “Y ( KEEN ) -S I GMAY ( KEEN ) 

YBU ( I > “ Y ( KEEN ) + S I GMAY ( KEEN ) 

YBU ( I + 1 > “Y ( KEEN ) +S I GMAY ( KEEN ) 

YBU( 1+2) “Y< KEEN) +S I GMAY (KEEN) 

YBU ( I +3 > =Y ( KEEN ) +SI GMAY ( KEEN ) 

END IF 

420 CONTINUE 
C 

IFUNSL.EQ. 1) CALL YTYPE<2) 

IF( INSL. EG. 2) CALL XTYPE(2) 

IF ( IANS. EG. 3) CALL YTYPE(2) 

IF( IANS. EG. 3) CALL XTYPE(2) 

CALL ZLINE(-4) 

CALL SYMBL ( ISYMB ) 

CALL XFRM(3) 

CALL XMFRM( 3) 

CALL YFRM(3) 

CALL YMFRM(3) 

IF(NSETX. EG. 1) CALL XNEAT(O) 

IF ( NSETY. EG. 1 ) CALL YNEAT(O) 

IF(NSETX. EG. 1) CALL DLIMX ( XMIN, XMAX ) 

IF (NSETY. EG. 1) CALL DLIMY(YMINj YMAX) 

CALL CHECMXDATA. YDATA) 

CALL DSPLAY(XDATA, YDATA) 

C 

IF ( MODE. EG. 1 ) THEN 

CAL L ERRBAR ( XDA TA, YBU/ YBD ) 

END IF 
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X1=X< 1 ) 

X2=X(NPTS> 

XINC=(X2-X1 > /FLOAT <4*NPTS> 

IN=0 

DO 430 XV=X1* X2, XINC 

IN=IN+1 

X< IN)=XV 

YF I T ( I N ) =FUNCTN ( X > IN. A) 

430 CONTINUE 

FILL DUMMY ARRAY WITH FITTING FUNCTION 

DO 440 1=2, IN 
XDATA( I )=X ( I ) 

YDATA< I )=YFIT< I ) 

440 CONTINUE 

CALL ZLINE(O) 

CALL SYMBL(O) 

CALL CPLOT < XDATA, YDATA) 

LABELING AXES 

CALL PRETTY (YLABEL, XLABEL. IYLEN, IXLEN, IXLAB, IYLAB ) 

IVY = IFIX< <575. -13. *IYLEN>/2. )+125 
CALL KA12AS< 50, YLABEL, IYLAB) 

CALL CHAR <20, IVY, IYLAB, 50, 90. , 1. ) 

IVX = IFIX< (750. -13. *IXLEN>/2. >+150 
CALL KA12AS( 50, XLABEL, IXLAB) 

CALL CHAR( IVX, 20, IXLAB, 50, 0. , 1. ) 

CALL FRAME 
CALL BELL 
CALL T INPUT ( I ) 

CALL ERASE 

CALL FINITT<0, 700) 

STOP 

240 CONTINUE 
TRACK=0. 0 
WRITE (6, 450) 

DO 460 1= 1 , NPTS 

DIFF=< <Y< I )-YFIT ( I ) )/Y< I ) )*100. 0 
TRACK=TRACK+DIFF 

WRITE <6, 470) X < I ) , Y ( I ) , YFIT < I ) , DIFF 
460 CONTINUE 

TRACK-TRACK/FLOAT (NPTS) 

WRITE (6, 480) TRACK 
GO TO 490 

10 FOR MAT ( / , 'NONLINEAR CURV E-F ITTING CODE',/) 
20 FOR MAT (13) 

40 FOR MAT ( /, 'TOO MANY DATA POINTS IN RAWDAT, CHECK PARAMETER',/) 

50 FOR MAT ( / , 'NUMBER OF DATA PAIRS =',I3> 

60 FOR MAT ( ! , 'CHOSEN FITTING FUNCTION IS: ',/) 
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70 FORMAT (/ 'ENTER INITIAL GUESSES FOR THE A1 — >A',U, ' PARAMETERS'./) 
90 FOR MAT < 'FOR A(',Ii. ') ENTER GUESS?') 

150 FORMAT( 'FINISHED ITERATION #',I2, ' WITH REDUCED CHI. SG. - ', 1PE13. 4) 
175 FORMAT ( 'ITERATION STOPPED BECAUSE ABS < XCHI ) . LT. 0. 01 ' ) 

260 FORMAT( 'INPUT TITLE OF Y - AXIS') 

280 FORMAT( 'INPUT TITLE OF X - AXIS') 

270 FORMAT ( 80A1 ) 

110 FORMAT <5X, 'A(', II, ')*=', 1PE14. 6) 

100 FORMAT </, 'STARTING VALUES',/) 

170 FORMAT(/, 'THERE WERE ',13,' ITERATIONS') 

190 FORMAT (/, 'USING ') 

45 FOR MAT ( 'THE FINAL COEFFICIENTS ARE',/) 

210 FORMAT < 5X» 'A(', II, ' )«* ', 1PE14. 6) 

220 FORMAT (/, 'WITH REDUCED CHI. SQUARE*®', 1PE16. 6) 

230 FORMAT*/, 'DO YOU WANT A DATA REVIEW ?<1«YES, 0«N0> ' ) 

250 FORMAT < /, 'DO YOU WANT TO PLOT DATA ?<1-YES, 0-N0>',/) 

480 FOR MAT ( / , 'MEAN OF 7. ERROR = ',F14. 6) 

450 FORMAT < / , 5X > 'X-DATA', 14X, 'Y-DATA ' , 13X, 

- 'YF IT ', 14X, '7. DIFFR. ') 

470 FORMAT* IX, 3<1PE14. 6, 5X), 1PE14. 6) 

290 FOR MAT < 'WHICH TYPE OF GRAPH DO YOU WANT? ', //, 5X, 

- '1 - LINEAR',/, 5X, '2 - SEMI-LOG ', /, 5X, '3 - LOG-LOG', 

- //, 'INPUT THE NUMBER OF YOUR SELECTION ?') 

320 FOR MAT < / , 'DO YOU WANT: ', //5X, '1 - LOG Y',/5X, '2 - LOG X',// 

- , 'INPUT WHICH <1=Y, 2“X>? ' ) 

330 FORMAT*/, 'DO YOU WANT SPECIAL SYMBOLS TO DENOTE DATA POINTS', 

- /, ' <1=YES» 0®*I0>?'> 

350 FOR MAT < 'SYMBOLS ARE: ', //, 

- 6X, '1 - CIRCLE', /, 

- 6X, '2 - CROSS', /, 

- 6X, '3 - TRIANGLE', /, 

- 6X, '4 - SQUARE',/, 

- 6X, '5 - STAR', /, 

- 6X, '6 - DIAMOND ', /, 

- 6X, '7 - VERTICAL BAR',/, 

- 6X, '8 - + SYMBOL', /, 

- 6X, '9 - UP ARROW BELOW POINT',/, 

- 5X, '10 - DOWN ARROW BELOW POINT',/, 

- 5X, '11 - REVERSE TRIANGLE',//, 

- 'INPUT THE NUMBER MATCHING YOUR SELECTION ?') 

360 FOR MAT ( / , 'DO YOU WANT TO SET THE X RANGE',/. '<1*®YES, 0-N0>? ' ) 

380 FORMAT (/, 'INPUT XMIN, XMAX ?') 

390 FORMAT (/, 'DO YOU WANT TO SET THE Y RANGE',/, 'Cl-YES, 0-N0>? ' ) 

410 FORMAT (/, 'INPUT YMIN, YMAX ?') 

120 FOR MAT < / ) 

END 
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c 

C********************** FUNCTION FUNCTN< X, I, A) ****************** 
PURPOSE 

EVALUATE TERMS OF FUNCTION FOR NON-LINEAR LEAST-SQUARES SEARCH 
USAGE 

RESULT=FUNCTN< X , I, A) 

DESCRIPTION OF PARAMETERS 

LL - NO. OF COEFFICIENTS OF FITTING FUNCTION 
X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE 
I - INDEX OF DATA POINT 
A - ARRAY OF PARAMETERS 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 

FUNCTION FUNCTNIX, I, A) 

PAR AMETER < LL*9 ) 

DIMENSION X(l), A(LL) 

XI-XU ) 

FUNCTN-AU )*TEXP(-A(2)*XI )*COS( A(3)*XI+A<4) >+A<5> 

FUNCTN*A< 1 )*TEXP<-0. 5*( ( X I-A<2) ) /A( 3) )**2) + 

A<4 ) *TEXP ( -0. 5*( (XI-A(5) >/A(6> )**2)+ 

A < 7 ) +A( 8)*XI+A<9)*XI *#2 

RETURN 
END 


***************** FUNCTION TEXP(X) *************************** 
PURPOSE 

TO ELIMINATE OVER/UNDER FLOW OF CPU IF EXP IS USED 
USAGE 

TEXP*EXP(X> 

FUNCTION TEXP(X) 

C 

IF ( X . LT. -100. ) X=-100. 

IF ( X . GT. 100. ) X=100. 

TEXP=EXP(X) 

END 
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c 

C**+* SUBROUTINE PRETTY( YLABEL, XLABEL, IYLEN, IXLEN, IXLAB, IYLAB ) ***** 
PURPOSE 

TO SEARCH THROUGH X AND Y TITLE AND COUNT THE NUMBER OF CHARACTERS 
USAGE 

CALL PRETTY (YLABEL, XLABEL, IYLEN. IXLEN, IXLAB. IYLAB) 

DESCRIPTION OF PARAMETERS 
YLABEL - ARRAY OF Y-TITLE 
XLABEL - ARRAY OF X-TITLE 
IYLEN - LENGTH OF ARRAY YLABEL 
IXLEN - LENGTH OF ARRAY XLABEL 

IYLAB - ARRAY CONTAINING INTEGER EQUIVALENT OF YLABEL 
IXLAB - ARRAY CONTAINING INTEGER EQUIVALENT OF XLABEL 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
CALL KA12AS( 50, YLABEL, IYLAB) : KA12AS IS A PLOT-10 CALL 

SUBROUTINE PRETTY ( YLABEL, XLABEL, IY. IX, IXLAB, IYLAB) 

C 

DIMENSION YLABEL (1), XLABEL < 1 ), IXLAB(l), IYLAB(l) 

C 

NCHAR=32 

CALL KA 1 2AS (■ 50, Y LABEL , I YLAB ) 

CALL KA 1 2AS ( 50, XLABEL, IXLAB) 

C 

DO 10 1 = 1, 50 

IF( ( IYLAB ( I ) . EQ. NCHAR ) . AND. (IYLAB(I + 1). EQ. NCHAR ) . AND. 

- ( IYLAB ( I +2 ) . EQ. NCHAR ) ) IY=I 
IF< ( IYLAB < I ) . EQ. NCHAR). AND. ( IYLAB ( 1 + 1 ) . EQ. NCHAR). AND. 

-< IYLAB ( 1+2). EQ. NCHAR)) GO TO 20 
10 CONTINUE 
C 

20 DO 30 1=1, 50 

IF< ( IXLAB(I). EQ. NCHAR). AND. (IXLAB<I + 1). EQ. NCHAR). AND. 

-(IXLAB (1+2). EQ. NCHAR)) IX=I 
IF( ( IXLAB(I). EQ. NCHAR). AND. (IXLAB(I + 1). EQ. NCHAR). AND. 

-( IXLAB( 1+2). EQ NCHAR) ) GO TO 40 
30 CONTINUE 
C 

40 RETURN 
END 
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c 

c ******* SUBROUTINE CHAR (LOCX, LOCY, ISTRNG, NCHAR, ANG, SIG) ********* 
PURPOSE 

REQUIRED PLOT-1 0 SUBROUTINE < CHARACTER GENERATION PACKAGE) 
USAGE 

CALL CHAR < LOCX< LOCY, ISTRNG, NCHAR, ANG, SIG) 

DESCRIPTION OF PARAMETERS 

LOCX - INTEGER VALUE OF LOCATION OF XDOT (0 — >1024) 

LOCY - INTEGER VALUE OF LOCATION OF YDOT <0~>780) 

ISTRNG - ARRAY CONTAINING TITLE STRING 
NCHAR - NO. OF CHARACTER IN ISTRNG 
ANG - ROTATION ANGLE FOR PLOTTING CHARACTER 
SIG - SIZE OF PLOTTED CHARACTER 

SUBROUTINE CHAR < LOCX, LOCY, ISTRNG, NCHAR, ANG, SIZ ) 

C 

DIMENSION ISTRNG <NCHAR) 

C 

REAL COMST < 60 ) 

C 

CALL SVSTAT( COMST) 

CALL RESET 
CALL BINITT 

CALL MO VABS< LOCX, LOCY) 

CALL RROTAT ( ANG ) 

CALL ZR5CALE(SIZ ) 

C 

DO 10 I“l» NCHAR 
CALL LCHAR < ISTRNG ( I ) ) 

10 CONTINUE 
C 

CALL RESTAT( COMST) 

RETURN 

END 
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c 

C********** SUBROUTINE ERRBAR < XDATA. YBU, YBD) *************** 
PURPOSE 

TO DRAW ERROR BAR IF M0DE=1 
USAGE 

CALL ERRBAR < XDATA. YBU. YBD ) 

DESCRIPTION OF PARAMETERS 

XDATA - DUMMY ARRAY TO STORE INDEPENDENT DATA POINTS 
YBU - DUMMY ARRAY TO STORE SIGMAY 
YBD - DUMMY ARRAY TO STORE SIGMAY 
SUBROUTINE ERRBAR < XDATA. YBU. YBD) 

DIMENSION XDATA ( 1 ), YBU<1). YBD<1) 

NDATA-XDATA(l) 

XMIN-1. E+99 
XMAX— XMIN 

DO 5 1-2, NDATA+ 1 

IF < XDATA ( I ) . LT. XMIN) XMIN-XDATA( I ) 

IF ( XDATA ( I ) . GT. XMAX > XMAX-XDATA < I ) 

CONTINUE 

WEERB- ( XMAX-XMI N > / < 2. 0*FL0AT ( NDATA ) ) 

C 

DO 10 1-2, NDATA+1 
XL- XDATA ( I ) -WEERB 
XR- XDATA ( I ) +WEERB 
CALL M0VEA< XL. YBU( I ) ) 

CALL DRAWA< XR» YBU< I ) ) 

10 CONTINUE 

C 

DO 15 1-2, NDATA + 1 

CALL MO VEA ( XDAT A ( I ) , YBU ( I ) ) 

CALL DR AWA < XDAT A ( I ) , YBD ( I ) ) 

15 CONTINUE 

C 

DO 20 1-2, NDATA + 1 
XL- XDATA < I > -WEERB 
XR- XDATA ( I ) +WEERB 
CALL MOVEA< XL, YBD( I ) ) 

CALL DR AWA ( XR. YBD< I ) ) 

20 CONTINUE 

C - 

RETURN 

END 
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c 

c**** SUBROUTINE CHIFI T< X, Y, SIGMAY, NPTS, NTERMS, MODE, A, DELTAA, SIGMAA, 

C - YFIT, CHISQR) ^************************* 

C 

C PURPOSE 

C MAKE A LEAST-SQUARES FIT TO A NON-LINEAR FUNCTION 

C WITH A PARABOLIC EXPANSION OF CHI. SQUARE 

C 

C SOURCE 

C DATA REDUCTION AND ERROR ANALYSIS FOR THE PHYSICAL SCIENCES 

C P. R. BEVINGTON 

C 

C USAGE 

C CALL CHIFIT ( X. Y» SIGMAY, NPTS. NTERMS, MODE, A, DELTAA, SIGMAA, YFIT. CHISQR ) 

C 

C DESCRIPTION OF PARAMETERS 

C LL - NO. OF COEFFICIENTS OF FITTING FUNCTION 

C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE 

C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE 

C SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS 

C NPTS - NUMBER OF PAIRS OF DATA POINTS 

C NTERMS - NUMBER OF PARAMETERS 

C MODE - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT 

C +1 (INSTRUMENTAL) WEIGHT < I ) = 1 . /SIGMAY ( I >**2 

C 0 (NO WEIGHTING) WEIGHT ( I )«1. 0 

C -1 (STATISTICAL) WE IGHT ( I ) = 1 . /Y ( I ) 

C A - ARRAY OF PARAMETERS 

C DELTAA - ARRAY OF INCREMENTS FOR PARAMETERS A 

C SIGMAA - ARRAY OF STANDARD DEVIATIONS FOR PARAMETERS A 

C YFIT - ARRAY OF CALCULATED VALUES OF Y 

C CHISQR - REDUCED CHI. SQUARE FOR FIT 

C 

C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 

C 

C FUNCTN(X, I, A) 

C EVALUATES THE FITTING FUNCTION FOR THE I-TH TERM 

C 

C FCHI SQ( Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 

C EVALUATES REDUCED CHI. SQUARE FOR FIT TO DATA 

C 

C MATINV(ARRAY, NTERMS, DET) 

C INVERTS A SYMMETRIC TWO-DIMENSIONAL MATRIX OF DEGREE NTERMS 

C AND CALCULATES ITS DETERMINANTS 

C 

C COMMENTS 

C DIMENSION STATEMENT VALID FOR NTERMS IS CHANGED BY PARAMETER STATEMENT 

C 

SUBROUTINE CHIFIT(X, Y, SIGMAY, NPTS, NTERMS, MODE, A, DELTAA, SIGMAA, 

-YFIT, CHISQR) 

C 

PAR AMETER ( LL"9 ) 

C 

DOUBLE PRECISION ALPHA 
C 

DIMENSION X < 1 ), Y ( 1 ) , SIGMAY ( 1 ) , A( 1 ) , DELTAA ( 1 ) , SIGMAA( 1 ) , YFIT ( 1 ) 
DIMENSION ALPHA ( LL, LL ) , BET A < LL ) , DA ( LL ) 
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11 NFREE-NPTS-NTERMS 
FREE -NFREE 
IF(NFREE) 14,14,16 
14 CHISGR-O. 

GO TO 120 

16 DO 17 1=1, NPTS 

17 YFIT ( I )=FUNCTN< X, I, A) 

18 CHISQ1-FCHISQ<Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 
EVALUATE ALPHA AND BETA MATRICES 

20 DO 60 J= 1 , NTERMB 
A(J) +DELTAA( J) 

21 AJ=A(J) 

A(U)=AJ+DELTAA<U) 

DO 24 1 = 1, NPTS 
24 YFIT ( I )=FUNCTN< X , I , A) 

CHISQ2=FCHISQ<Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 
ALPHA (J, J)=CHISG2-2. *CHISQ1 
BETA < J)=-CHISG2 

31 DO 50 K=1 , NTERMS 
IF(K-J) 33,50,36 

33 ALPHA (K, J) = <ALPHA(K, J) -CHISQ2 ) /2. 

ALPH A ( J, K ) -ALPHA ( K, J ) 

GO TO 50 

36 ALPHA< J, K)=CHISG1-CHISQ2 

A ( J ) +DELT AA ( J ) AND A(K)+DELTAA(K) 

41 AK=A ( K ) 

A ( K ) “AK+DELT AA ( K ) 

DO 44 1 = 1, NPTS 
44 YFIT < I )=FUNCTN< X, I, A) 

CHISQ3=FCHISQ(Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 

ALPHAt J, K ) =ALPHA ( J, K)+CHISQ3 

A(K)=AK 

50 CONTINUE 

A< J) -DELTAA( J) 

51 A( J) =AU-DELTAA< V ) 
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DO 53 I-l > NPTS 
53 YFI T ( I ) =FUNC TN< X , I » A ) 

CHISQ3=FCHISQ<Y, SIGMAY, NPTS, NFREE, MODE, YFIT ) 
A<J)=AJ 

AI_PHA( J, J) = < ALPHA< J, J)+CHISQ3> /2. 

BETA < J)= (BETA< J) +CHISQ3) /4. 

60 CONTINUE 

ELIMINATE NEGATIVE CURVATURE 

61 DO 70 J= 1 , NTERMS 
IF<ALPHA<U< J) ) 63,65,70 

63 ALPHA< J, J)=-ALPHA< J, J) 

GO TO 66 

65 ALPHA < J, J)=0. 01 

66 DO 72 K-l, NTERMS 
IF<K-U) 68,72,68 

68 ALPHA <U, K)=0. 0 
ALPHA <K, J)»0. 0 
72 CONTINUE 

70 CONTINUE 

INVERT MATRIX AND EVALUATE PARAMETER INCREMENTS 

71 CALL MATINV<ALPHA, NTERMS, DET) 

DO 76 J-l, NTERMS 
DA<J)-0. 0 

74 DO 75 K-l, NTERMS 

75 DA< J )«DA < U)+BETA < K) *ALPHA< J, K ) 

76 DA<U )=0. 2*DA < J) *DELTAA < J ) 

MAKE SURE CHI. SQUARE DECREASES 

81 DO 82 J-l, NTERMS 

82 A<J)-A<J)+DA< J) 

C 

83 DO 84 1-1, NPTS 

84 YFIT ( I )-FUNCTN<X, I, A) 

C 

CHISQ2«FCHISQ<Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 
IF<CHISQ1-CHI 9Q2 ) 87,91,91 
C 

87 DO 89 J-l, NTERMS 
DA< J )=DA < J)/2. 

89 A ( J) -A< J >~DA< J) 
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c 


GO TO S3 


INCREMENT PARAMETERS UNTIL CHI. SQUARE STARTS TO INCREASE 

91 DO 92 J= 1 , NTERMS 

92 A<J)=A< J)+DA<J) 

DO 94 1 = 1. NPTS 
94 YFIT ( I )=FUNCTN< X. I. A) 

CHISQ3=FCHISQ(Y. SIGMAY. NPTS. NFREE, MODE. YFIT) 

IF(CHISQ3— CHISQ2 ) 97.101,101 
97 CHISQ1=CHISQ2 
CHISQ2-CHISQ3 
99 GO TO 91 

FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS 

101 DELTA* 1 1. /(I. -*-<CHISQl— CHISQ2)/ (CHISG3— CHISQ2) )+. 5 

DO 104 J=l. NTERMS 
A(J) =A< J)-DELTA*DA( J) 

104 S IGM AA ( J ) =DELTAA ( J ) *SGRT ( FREE*DABS ( ALPHA < J. J) ) ) 

DO 106 1 = 1, NPTS 
106 YFIT < I )=FUNCTN< X, I, A) 

CHISQR=FCHISQ<Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 

111 IF<CHISQ2— CHI9QR ) 112,120,120 
C 

112 DO 113 J=l. NTERMS 

113 A(J)=A(J) + (DELTA-1. )*DA<J) 

C 

DO 115 1 = 1, NPTS 
115 YFIT ( I )=FUNCTN< X , I, A) 

C 

CHISQR=CHISQ2 
120 RETURN 
END 
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c 

C*« ******** FUNCTION FCHISQ(Y, SIGMAY, NPTS, NFREE, MODE, YFIT) ********* 
PURPOSE 

EVALUATE REDUCED CHI. SQUARE FOR FIT TO DATA 
FCH I SQ-SUM < ( Y-Y F I T ) **2/S I GMA**2 ) /NFREE 

SOURCE 

DATA REDUCTION AND ERROR ANALYSIS FOR THE PHYSICAL SCIENCES 
P. R. BEVINGTON 

USAGE 

RESULT=FCHISQ(Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 

DESCRIPTION OF PARAMETERS 
Y - ARRAY OF DATA POINTS 

SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR DATA POINTS 
NPTS - NUMBER OF DATA POINTS 
NFREE - NUMBER OF DEGREES OF FREEDOM 

MODE - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT 
+ 1 (INSTRUMENTAL) WEIGHT< I ) = 1. /SIGMAY< I )**2 
0 (NO WEIGHTING) WEIGHT ( I ) = 1. 0 
-1 (STATISTICAL) WEIGHT( I ) = 1. /Y( I ) 

YFIT - ARRAY OF CALCULATED VALUES OF Y 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 

FUNCTION FCHISQ ( Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 

DOUBLE PRECISION CH I SQ, WEIGHT 
DIMENSION Y ( 1 ) , S I GMAY ( 1 ) , YF I T ( 1 ) 

11 CHISQ-O. 

12 IF( NFREE) 13,13,20 

13 FCH ISG=0. 0 
GO TO 40 

ACCUMULATE CHI. SQUARES 

20 DO 30 1 = 1, NPTS 

21 IF (MODE) 22,27,29 

22 IF ( Y ( I ) ) 25,27 , 23 

23 WEIGHT-1. /Y( I) 

GO TO 30 

25 WE I GHT = 1 . / ( - Y< I ) ) 

GO TO 30 
27 WEIGHT-1. 

GO TO 30 

29 WE I GHT = 1 . /SI GMA Y ( I ) **2 

30 CHI SQ-CHISQ-+WEI GHT* ( Y ( I )-YFIT ( I ) )**2 

DIVIDE BY NUMBER OF DEGREES OF FREEDOM 

31 FREE-NFREE 

32 FCH ISQ-CHISQ/FREE 
40 RETURN 

END 
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c 

c ************* SUBROUTINE MATINV ( ARRAY, NORDER, DET ) ************** 

C 

C PURPOSE 

C INVERT A SYMMETRIC MATRIX AND CALCULATE ITS DETERMINANT 

C 

C SOURCE 

C DATA REDUCTION AND ERROR ANALYSIS FOR THE PHYSICAL SCIENCES 

C P. R. BEVINGTON 

C 

C USAGE 

C CALL MATINV(ARRAY, NORDER, DET) 

C 

C DESCRIPTION OF PARAMETERS 

C LL - NO. OF COEFFICIENTS OF FITTING FUNCTION 

C ARRAY - INPUT MATRIX WHICH IS REPLACED BY ITS INVERSE 

C NORDER - DEGREE OF MATRIX (ORDER OF DETERMINANT) 

C DET - DETERMINANT OF INPUT MATRIX 

C 

C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 

C NONE 

C 

C COMMENTS 

C DIMENSION STATEMENT VALID FOR NORDER IS CHANGED BY PARAMETER STATEMENT 

C 

SUBROUTINE MAT I NV< ARRAY, NORDER, DET) 

C 

PAR AMETER ( LL=9 ) 

C 

DOUBLE PRECISION ARRAY, AMAX, SAVE 
C 

DIMENSION ARRAY (LL, LL), IK(LL). JK(LL) 

C 

10 DET-1.0 
C 

11 DO 100 K-l. NORDER 

FIND LARGEST ELEMENT ARRAY ( I , J) IN REST OF MATRIX 
AMAX«=0. 

21 DO 30 I =K, NORDER 

DO 30 J-K, NORDER 

23 IF ( DABS ( AMAX )-D ABS < ARRAY ( I , J ) ) ) 24,24,30 

24 AMAX=ARRAY( I, J) 

IK( K)*I 
JK< K)=J 

30 CONTINUE 

INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY (K, K) 

31 IF( AMAX) 41, 32, 41 

32 DET =0. 0 
GO TO 140 

41 1 = 1 K < K ) 

IF( I— K ) 21, 51, 43 
C 
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43 DO 50 J=l, NORDER 
SAVE=ARRAY<K» J) 

ARRAY<K» J)=ARRAY< I, J) 

50 ARRAY (I, J)=-SAVE 
C 

51 J=JK<K> 

IF(J-K) 21,61,53 

C 

53 DO 60 1 = 1, NORDER 
SAVE=ARRAY ( I , K ) 

ARRAY ( I , K)=ARRAY < I, J) 

60 ARRAYd, J)=-SAVE 

ACCUMULATE ELEMENTS OF INVERSE MATRIX 

61 DO 70 1 = 1, NORDER 
IF(I-K) 63,70,63 

63 ARRAYd, K)=-ARRAY< I, K)/AMAX 

70 CONTINUE 

71 DO 80 1 = 1, NORDER 
DO 80 J=l, NORDER 
IF < I— K ) 74,80,74 

74 IF(J-K) 75,80,75 

75 ARRAYd, J)=ARRAYd, J)+ARRAY<I, K)*ARRAY(K, J) 

80 CONTINUE 

I 

81 DO 90 J=l, NORDER 
IF(J-K) 83,90,83 

83 ARRAY<K. J)=ARRAY<K, Jl/AMAX 
90 CONTINUE 

ARRAY<K, K) = l. /AMAX 

100 DET =DET *AMAX 

RESTORE ORDERING OF MATRIX 

101 DO 130 L=l, NORDER 
K=N ORDER -L+l 
J=IK<K> 

IF(J-K) 111,111,105 

105 DO 110 1=1, NORDER 
SAVE=ARRAY(I,K) 

ARRAYd, K)=-ARRAY(I, J) 

110 ARRAYd, U)=SAVE 

111 I=JK(K) 

IF(I-K) 130.130,113 

113 DO 120 J*l, NORDER 
SAVE=ARRAY < K, J > 

ARRAY <K, J)=-ARRAY(I, J) 

120 ARRAYd, J)=SAVE 
C 

130 CONTINUE 
C 

140 RETURN 
END 
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Appendix B 

Sample Case 1: Application in Classical Mechanics 

Appendix B is the complete listing of an interactive session for a fitting function of five parameters. 
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NONLINEAR C U R V E-F ITTING CODE 


NUMBER OF DATA PAIRS =166 
CHOSEN FITTING FUNCTION IS: 

Y = A*EXP(-B*X)*COS(C*X+D)+E 

ENTER INITIAL GUESSES FOR THE Al — >A5 PARAMETERS 


FOR A ( 1 ) 
? 68.0 
FOR A( 2 ) 
? 0.002 
FOR A( 3) 
? 0.1 
FOR A( 4 ) 
? - 2.0 
FOR A( 5) 
? 100.0 


ENTER GUESS? 
ENTER GUESS? 
ENTER GUESS? 
ENTER GUESS? 
ENTER GUESS? 


STARTING VALUES 


AC 1 ) = 6 . 300000E+01 

A( 2) = 2 . 0000 OOE-O 3 

A( 3) = 1 .OOOOOOE-Ol 

AC4>= -2 . OOOOOQE+OO 
AC5)= 1 .OOOOOOE+02 


FINISHED ITERATION # 1 WITH REDUCED CHI.SQ.= 2 

FINISHED ITERATION # 2 WITH REDUCED CHI.SQ.* 1 

FINISHED ITERATION # 3 WITH REDUCED CHI. SO. = 4 

FINISHED ITERATION # 4 WITH REDUCED CHI.SQ.= 2 

FINISHED ITERATION # 5 WITH REDUCED CHI.SQ.= 1 

FINISHED ITERATION # 6 WITH REDUCED CHI.SQ.= 1 

ITERATION STOPPED BECAUSE ABSCXCHI ) . LT . 0 . 01 

THERE WERE 6 ITERATIONS 

USING 

Y = A*EXPC~B*X)*COSCC*X+D)+E 
THE FINAL COEFFICIENTS ARE 

A( 1 ) = 7 . 3776S4E+01 

AC 2) = 1 .140380E-03 

AC 3) = 9.244445E-02 

AC 4) = -3 . 744320 E+GO 
AC 5) = 1 . 343607E+02 

WITH REDUCED CHI. SQUARE= 1.232234E-01 


3147E+01 
3383E+01 
. 363QE4-Q0 
0090E-01 
. 2535E-Q 1 
2322E-01 
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DO YOU WANT A DATA REVIEW ?<1=YES, 0=NO> 

? 0 

DO YOU WANT TO PLOT DATA ?<1=YES, 0=N0> 

? 1 

INPUT TITLE OF Y - AXIS 
amplitude (cm.) 

INPUT TITLE OF X - AXIS 
? time (sec.) 

WHICH TYPE OF GRAPH DO YOU WANT? 

1 - LINEAR 

2 - SEMI -LOG 

3 - LOG-LOG 

INPUT THE NUMBER OF YOUR SELECTION ? 

? 1 

DO YOU WANT SPECIAL SYMBOLS TO DENOTE DATA POINTS 
<1=YES, G=NO> ? 

? 1 

SYMBOLS ARE: 

1 - CIRCLE 

2 - CROSS 

3 - TRIANGLE 

4 - SQUARE 

5 - STAR 

6 - DIAMOND 

7 - VERTICAL BAR 
3 - + SYMBOL 

9 - UP ARROW BELOW POINT 

10 - DOWN ARROW BELOW POINT 

11 - REVERSE TRIANGLE 

INPUT THE NUMBER MATCHING YOUR SELECTION ? 

? 1 

DO YOU WANT TO SET THE X RANGE 
<1=YES , 0=NO>? 

* 0 

DO YOU WANT TO SET THE Y RANGE 
<1 =YES , 0=NO>? 

? 1 

INPUT YMIN, YMAX ? 

? 0.0,250.0 


Appendix C 

Sample Case 2: Application in Fluid Mechanics 

Appendix C is the complete listing of an interactive session for a fitting function of nine parameters. 
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nonlinear guru e-f itting code 


NUMBER OF DATA PAIRS = 22 
CHOSEN FITTING FUNCTION IS: 

y = A*EXPC -0 .5*C CX-B)/C)**2)+D*EXP( -0 . 5*C (X-E)/F )^2)+G+H*X+I*X**2 

ENTER INITIAL GUESSES FOR THE Al — >A9 PARAMETERS 

FOR Ad) ENTER GUESS? 

? - 1.0 

FOR A( 2) ENTER GUESS? 

? 0.0 

FOR A( 3) ENTER GUESS? 

? 3.5 

FOR AC 4 ) ENTER GUESS? 

? 0 .5 

FOR A( 5) ENTER GUESS? 

? 6.0 

FOR AOS) ENTER GUESS? 

? 3.0 

FOR A( 7) ENTER GUESS? 

? -1 .0 

FOR ACS) ENTER GUESS? 

? -0.5 

FOR ACS) ENTER GUESS? 

? - 0.1 

STARTING VALUES 

A C 1 ) =’ -1 .000000E+00 
A C 2 ) = . 000000E+00 

AC 3) = 3 . 500 00 0E+0 0 

A ( 4 ) = 5 . OOOOOOE-Ol 

AC5)= 6 . 000000E+00 

AC 6) = 3 . OGOGQOE+GQ 

AC 7 ) = -1 .000000E+00 
AC 8) = -5. OOOOOOE-Ol 
A ( 3 ) = -1 .OOOOOOE-Ol 


FINISHED ITERATION # 1 WITH REDUCED CHI.SQ.= 3.3716E+03 
FINISHED ITERATION # 2 WITH REDUCED CHI. SO. = 1.6085E+00 
FINISHED ITERATION # 3 WITH REDUCED CHI. SO. = 7.0243E-01 
FINISHED ITERATION # 4 WITH REDUCED CHI.SQ.= 6.9253E-01 
ITERATION STOPPED BECAUSE ABSCXCHI ) .LT.O .01 
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THERE WERE 4 ITERATIONS 


USING 

y = A*EXP(-0 .5*(<X-B)/C)**2)+D*EXP(-0 .5*( (X-E)/F)**2)+G+H*X+I*X**2 
THE FINAL COEFFICIENTS ARE 

A ( 1 ) = -6 • 533788E-0 2 
A( 2) = .OQGOOOE+OO 
A ( 3 ) = 3.431341E+00 

A( 4) = 2.952222E-02 

A ( 5 ) = 5 • 160242E+00 

A(S> = 2 . 575265E+0O 

A ( 7 ) = -2 . 212542E-02 
A ( 8 ) = -4 . 082130E-03 
A ( 9 ) = -2 . 473477E-G5 


WITH REDUCED CHI. SQUARE= 6.925316E-01 


DO YOU WANT A DATA 

9 i 

REVIEW ?<1=YES , 

O 

Z 

ll 

o 


X-DATA 

Y-DATA 

YFIT 

% DIFFR. 

-3.00 0000 E+0 1 

6 . 7660 00 E-0 2 

7 . 80771 8E-02 

-1 .539637E+01 

-2.2019G0E+01 

5.803000E-02 

5.576663E-02 

3 . 932988E+00 

-1 .901710E+Q1 

5.212000E-02 

4 . 855950 E-0 2 

1 .066864E+01 

-1 .601690E+Q1 

4.021000E-02 

3 . 691 094E-02 

3 . 204533E+GG 

-1 .301570E+01 

2 . 878000 E-0 2 

2.676663E— 02 

6. 99571 5E+00 

-1 . 00 1340E+0 1 

1 .427000E-02 

1 .533332E-02 

-7 . 489994E+00 

-S.01530GE+00 

5 . 51000GE-03 

4 . 706874E-0 3 

1 .457530E+01 

-6 . 0 0 570 0 E+G 0 

-1 .124000E-02 

-1 .273132E-02 

-1 . 326794E+Q 1 

-5 . 0 G620GE+Q 0 

-2.505000E-02 

-2 . 500952E-02 

1 .616062E-01 

-4 . OQ59GOE+00 

-4.011000E-02 

-3 . 94231 7E-02 

1 .71236QE+00 

-3.006100E+00 

-5.841000E-02 

-5 . 473962E-0 2 

6 . 28331 5E+00 

-2.005300E+00 

-7 . 1 23000 E -02 

-6.392002E-02 

3.242934E+00 

-1 .G060Q0E+0Q 

-8 . G43GG0E-Q2 

-7 . 943255E-0 2 

1 . 30150 4E+00 

-7 . 3000 00E-03 

-3.273000E-02 

-3 . 393985E-0 2 

-1 .522850E+00 

9.941G00E-01 

-7 . 949000 E-0 2 

-8.136233E-02 

-2.356Q62E+00 

1 .99540 0E+00 

-7 . 122000E-0 2 

-7 . 20 91 65E-0 2 

-1 .223338E+00 

2 . 99480 0E-R00 

-6 . 2220 00 E-0 2 

-5 . 382644E-02 

5 • 454124E+GQ 

4.9952G0E+Q0 

-4 . 1 350 00 E-0 2 

-3 . 6490 70 E-0 2 

1 .175163E+01 

6 . 99540 OE+OO 

-3 . 2890 00 E-0 2 

-3 . 723041E-02 

-1 .319675E+G1 

3.99670GE+00 

-4 . 212000E-02 

-5 . 323740E-02 

-2.639458E+01 

1 .09S690E+01 

-7 . 369000E-02 

-6 . 806898E-02 

7.627933E+00 

1 . 29957GE+0 1 

-7 . 090 000 E-0 2 

-7 .91151 5E-02 

-1 .158695E+Q1 

MEAN OF % ERROR = 

-.475987 
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DO YOU WANT TO PLOT DATA ?<1=YES, 0=NO> 

? 1 

INPUT TITLE OF Y - AXIS 
- pressure coef . 

INPUT TITLE OF X - AXIS 
? di stance (in.) 

WHICH TYPE OF GRAPH DO YOU WANT? 

1 - LINEAR 

2 - SEMI -LOG 

3 - LOG-LOG 

INPUT THE NUMBER OF YOUR SELECTION ? 

? 1 

DO YOU WANT SPECIAL SYMBOLS TO DENOTE DATA POINTS 
< 1 =YES , 0 =N 0 > ? 

? 1 

SYMBOLS ARE: 

1 - CIRCLE 

2 - CROSS 

3 - TRIANGLE 

4 - SQUARE 

5 - STAR 

6 - DIAMOND 

? - VERTICAL BAR 

S - + SYMBOL 

3 - UP ARROW BELOW F 0 1 NT 

10 - DOWN ARROW BELOW POINT 

11 - REVERSE TRIANGLE 

INPUT THE NUMBER MATCHING YOUR SELECTION ? 

? 1 

DO YOU WANT TO SET THE X RANGE 
<1=YES , 0=NO>? 

? 0 

DO YOU WANT TO SET THE Y RANGE 
< 1 =YES , 0 =NQ > ? 
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